

#include "Elem_Solid_3D_HEXA8.h"
#include "GeomDataLagrange.h"
#include "SolutionDataSolid.h"
#include "QuadratureUtil.h"
#include "TimeFunction.h"
#include "stressrecovery.h"
#include "MyTime.h"
#include "TimeFunction.h"
#include "utilitiesmaterial.h"


using namespace std;


extern  std::vector<unique_ptr<TimeFunction> > timeFunction;
extern  MyTime                 myTime;
extern  bool  debug;



Elem_Solid_3D_HEXA8::Elem_Solid_3D_HEXA8()
{
  if (debug) cout << " constructor Elem_Solid_3D_HEXA8 \n\n";
  
  ndim   = 3;
  ndof   = 3;
  npElem = 8;
  nsize  = npElem*ndof;

  ELEM_SHAPE = ELEM_SHAPE_HEXA;
}




Elem_Solid_3D_HEXA8::~Elem_Solid_3D_HEXA8()
{
}





void Elem_Solid_3D_HEXA8::prepareElemData()
{
    ElementBase::prepareElemData();

    if(ELEM_FORMULATION == "Bbar")
    {
      Nc.resize(npElem);
      dNc_dx.resize(npElem);
      dNc_dy.resize(npElem);
      dNc_dz.resize(npElem);
    }

    return;
}




double Elem_Solid_3D_HEXA8::computeVolume(bool init)
{
    int  ii, jj, gp;
    double  dvol, Jac, param[3];

    VectorXd  N(npElem), dN_dx(npElem), dN_dy(npElem), dN_dz(npElem);

    double xNode[npElem], yNode[npElem], zNode[npElem], xx, yy, zz;
  
    for(ii=0;ii<npElem;ii++)
    {
      xNode[ii] = GeomData->NodePosOrigLocal[nodeInLocalElem[ii]][0];
      yNode[ii] = GeomData->NodePosOrigLocal[nodeInLocalElem[ii]][1];
      zNode[ii] = GeomData->NodePosOrigLocal[nodeInLocalElem[ii]][2];

      // xNode[ii] = GeomData->NodePosOrig[nodeNums[ii]][0];
      // yNode[ii] = GeomData->NodePosOrig[nodeNums[ii]][1];
      // zNode[ii] = GeomData->NodePosOrig[nodeNums[ii]][2];

    }

    vector<double>  gausspoints1, gausspoints2, gausspoints3, gaussweights;
    getGaussPointsHexa(nGP, gausspoints1, gausspoints2, gausspoints3, gaussweights);

    elemVol = 0.0;
    for(gp=0; gp<nGP; gp++)
    {
        param[0] = gausspoints1[gp];
        param[1] = gausspoints2[gp];
        param[2] = gausspoints3[gp];

        if(init)
          GeomData->computeBasisFunctions3D(CONFIG_ORIGINAL, ELEM_SHAPE_HEXA, param, nodeInLocalElem, &N(0), &dN_dx(0), &dN_dy(0), &dN_dz(0), Jac);
        else
          GeomData->computeBasisFunctions3D(CONFIG_DEFORMED, ELEM_SHAPE_HEXA, param, nodeInLocalElem, &N(0), &dN_dx(0), &dN_dy(0), &dN_dz(0), Jac);

         dvol = gaussweights[gp]*Jac;

         elemVol += dvol;
    }

    return  elemVol;
}



int Elem_Solid_3D_HEXA8::calcMassMatrix(MatrixXd& Mlocal, bool MassLumping)
{
    if( (Mlocal.rows() != nsize) || (Mlocal.cols() != nsize) )
      Mlocal.resize(nsize, nsize);
    Mlocal.setZero();

    // mass lumping - row-wise sum
    if(MassLumping)
    {
      double fact = MatlData->getDensity()*computeVolume(true)/npElem;

      for(int ii=0; ii<nsize; ii++)
      {
        Mlocal(ii,ii) = fact;
      }
    }
    else
    {
      int  nGPt, ii, jj, gp, TI, TIp1, TIp2, TJ, TJp1, TJp2;
      double  fact, dvol0, Jac, bb1, cc1, param[3];

      double rho0 = MatlData->getDensity();

      VectorXd  N(npElem), dN_dx(npElem), dN_dy(npElem), dN_dz(npElem);

      vector<double>  gausspoints1, gausspoints2, gausspoints3, gaussweights;
      if(ELEM_SHAPE == ELEM_SHAPE_TETRA)
      {
        nGPt=10;
        getGaussPointsTetra(nGPt, gausspoints1, gausspoints2, gausspoints3, gaussweights);
      }
      else if(ELEM_SHAPE == ELEM_SHAPE_WEDGE)
      {
        nGPt=9;
        getGaussPointsWedge(nGPt, gausspoints1, gausspoints2, gausspoints3, gaussweights);
      }
      else if(ELEM_SHAPE == ELEM_SHAPE_HEXA)
      {
        nGPt=27;
        getGaussPointsHexa(nGPt, gausspoints1, gausspoints2, gausspoints3, gaussweights);
      }
      else
      {
        throw runtime_error("Element type not defined...");
      }

      bool   axsy  = (GeomData->getModelBehaviour() == 3);

      for(gp=0; gp<nGPt; gp++)
      {
          param[0] = gausspoints1[gp];
          param[1] = gausspoints2[gp];
          param[2] = gausspoints3[gp];

          GeomData->computeBasisFunctions3D(CONFIG_ORIGINAL, ELEM_SHAPE, param, nodeInLocalElem, &N(0), &dN_dx(0), &dN_dy(0), &dN_dz(0), Jac);

          dvol0 = gaussweights[gp]*Jac;

          for(ii=0;ii<npElem;ii++)
          {
            bb1 = (dvol0*rho0)*N[ii];

            TI   = 3*ii;
            TIp1 = TI+1;
            TIp2 = TI+2;

            for(jj=0; jj<npElem; jj++)
            {
              TJ   = 3*jj;
              TJp1 = TJ+1;
              TJp2 = TJ+2;

              fact  = bb1*N[jj];

              Mlocal(TI,   TJ)    += fact ;
              Mlocal(TIp1, TJp1)  += fact ;
              Mlocal(TIp2, TJp2)  += fact ;
            }
          }
      } //gp
      //cout << " elemVol = " << elemVol << endl;
      //printMatrix(Klocal);  printf("\n\n\n");
    }

  return 0;
}


double  Elem_Solid_3D_HEXA8::calcCriticalTimeStep(bool flag)
{
  double  dtCric;

  return  dtCric;
}




int Elem_Solid_3D_HEXA8::calcStiffnessAndResidual(MatrixXd& Klocal, VectorXd& Flocal)
{
    if(ELEM_FORMULATION == "Displacement")
      ElementBase::calcStiffnessAndResidual_Formulation_Displacement_3D(Klocal, Flocal);
    else if(ELEM_FORMULATION == "Bbar")
      calcStiffnessAndResidual_Formulation_Bbar(Klocal, Flocal);
    else if(ELEM_FORMULATION == "Mixed")
      calcStiffnessAndResidual_Formulation_Mixed(Klocal, Flocal);

    return 0;
}






int Elem_Solid_3D_HEXA8::calcStiffnessAndResidual_Formulation_Bbar(MatrixXd& Klocal, VectorXd& Flocal)
{
    int   err, index, ii, jj, kk, gp, TI, TIp1, TIp2, TJ, TJp1, TJp2;

    VectorXd  N(npElem), dN_dx(npElem), dN_dy(npElem), dN_dz(npElem);

    double  detF, detFc, trFc, fact, dvol, dvol0, Jac, volstrain;
    double  bb1, bb2, bb3, bb4, bb5, cc1, cc2, cc3, cc4, cc5, fact1, fact2, fact3;
    double  param[3], bforce[3]={0.0,0.0,0.0}, force[3];
    double  veloCur[3], acceCur[3], sig[3];

    int    sss   = GeomData->getModelBehaviour();
    bool  finite = MatlData->isFiniteStrain();

    double rho0  = MatlData->getDensity() ;
    double rho   = rho0;

    //bforce[0]   = elmDat[6]*timeFunction[0]->getFactor() ;
    //bforce[1]   = elmDat[7]*timeFunction[0]->getFactor() ;
    //bforce[2]   = elmDat[8]*timeFunction[0]->getFactor() ;

    double af   = SolnData->td(2);
    double acceFact1 = SolnData->td(5);
    double acceFact2 = acceFact1;
    double FiniteFact = (finite == true) ? 1.0 : 0.0;


    double xNode[npElem], yNode[npElem], zNode[npElem], xx, yy, zz;
    for(ii=0;ii<npElem;ii++)
    {
      xNode[ii] = GeomData->NodePosOrigLocal[nodeInLocalElem[ii]][0];
      yNode[ii] = GeomData->NodePosOrigLocal[nodeInLocalElem[ii]][1];
      zNode[ii] = GeomData->NodePosOrigLocal[nodeInLocalElem[ii]][2];
      // xNode[ii] = GeomData->NodePosOrig[nodeNums[ii]][0];
      // yNode[ii] = GeomData->NodePosOrig[nodeNums[ii]][1];
      // zNode[ii] = GeomData->NodePosOrig[nodeNums[ii]][2];
    }

    MatrixXd  Bbar1(3,9), Bbar2(9,3), Ktemp(3,3), bc(3,9);
    MatrixXd  Cmat(9,9), Gc(3,9), F(3,3), Fc(3,3);
    Cmat.setZero();
    Gc.setZero();
    F.setZero();
    Fc.setZero();
    VectorXd  stre(9), Gxs(3);
    stre.setZero();


    int nGP=1;
    vector<double>  gausspoints1, gausspoints2, gausspoints3, gaussweights;

    // Calculate values at the centroid
    if(npElem == 6)
    {
      getGaussPointsWedge(nGP, gausspoints1, gausspoints2, gausspoints3, gaussweights);
    }
    else if(npElem == 8)
    {
      getGaussPointsHexa(nGP, gausspoints1, gausspoints2, gausspoints3, gaussweights);
    }
    else if(npElem == 10)
    {
      getGaussPointsTetra(nGP, gausspoints1, gausspoints2, gausspoints3, gaussweights);
    }
    else if(npElem == 18)
    {
      getGaussPointsWedge(nGP, gausspoints1, gausspoints2, gausspoints3, gaussweights);
    }
    else if(npElem == 27)
    {
      getGaussPointsHexa(nGP, gausspoints1, gausspoints2, gausspoints3, gaussweights);
    }


    param[0] = gausspoints1[0];
    param[1] = gausspoints2[0];
    param[2] = gausspoints3[0];

    GeomData->computeBasisFunctions3D(CONFIG_ORIGINAL, ELEM_SHAPE, param, nodeInLocalElem, &N(0), &dNc_dx(0), &dNc_dy(0), &dNc_dz(0), Jac);

    computeDefGradCur(dNc_dx, dNc_dy, dNc_dz, Fc);

    trFc  = Fc.trace();
    detFc = Fc.determinant();


    if(npElem == 6)
    {
      nGP = 2;
      getGaussPointsWedge(nGP, gausspoints1, gausspoints2, gausspoints3, gaussweights);
    }
    else if(npElem == 8)
    {
      nGP = 8;
      getGaussPointsHexa(nGP, gausspoints1, gausspoints2, gausspoints3, gaussweights);
    }
    else if(npElem == 10)
    {
      nGP = 4;
      getGaussPointsTetra(nGP, gausspoints1, gausspoints2, gausspoints3, gaussweights);
    }
    else if(npElem == 18)
    {
      nGP = 9;
      getGaussPointsWedge(nGP, gausspoints1, gausspoints2, gausspoints3, gaussweights);
    }
    else if(npElem == 27)
    {
      nGP = 27;
      getGaussPointsHexa(nGP, gausspoints1, gausspoints2, gausspoints3, gaussweights);
    }



    if(Klocal.rows() != nsize)
    {
      Klocal.resize(nsize, nsize);
      Flocal.resize(nsize);
    }
    Klocal.setZero();
    Flocal.setZero();

    for(gp=0; gp<nGP; gp++)
    {
        param[0] = gausspoints1[gp];
        param[1] = gausspoints2[gp];
        param[2] = gausspoints3[gp];

        GeomData->computeBasisFunctions3D(CONFIG_ORIGINAL, ELEM_SHAPE, param, nodeInLocalElem, &N(0), &dN_dx(0), &dN_dy(0), &dN_dz(0), Jac);

        dvol0 = gaussweights[gp]*Jac;
        dvol  = dvol0;

        xx = yy = zz = 0.0;
        for(ii=0;ii<npElem;ii++)
        {
          xx += N[ii]*xNode[ii];
          yy += N[ii]*yNode[ii];
          zz += N[ii]*zNode[ii];
        }

        computeDefGradCur(dN_dx, dN_dy, dN_dz, F);

        volstrain = F.trace() - 3.0;
        detF = F.determinant();

        if(detF < 0.0)
        {
          throw runtime_error("Negative Jacobian in the element");
        }

        //acceCur[0] = computeValueDotDotCur(0, N);
        //acceCur[1] = computeValueDotDotCur(1, N);
        //acceCur[2] = computeValueDotDotCur(2, N);

        if(finite)
        {
          GeomData->computeBasisFunctions3D(CONFIG_DEFORMED, ELEM_SHAPE, param, nodeInLocalElem, &N(0), &dN_dx(0), &dN_dy(0), &dN_dz(0), Jac);
          dvol = gaussweights[gp]*Jac;
        }

        // ADJUST F
        fact = (trFc - F.trace())/3.0;

        F(0,0) += fact ;
        F(1,1) += fact ;
        F(2,2) += fact ;

        MatlData->computeStressAndTangent(true, sss, F, F, pres, stre, Cmat, ivar, gp, myTime.dt);
        //if(err !=0)    return 1;
        //printMatrix(Cmat);


        // Calculate Stiffness and Residual
        //==============================================

        force[0] = rho0*bforce[0];
        force[1] = rho0*bforce[1];
        force[2] = rho0*bforce[2];

        force[0] -= rho0*acceCur[0];
        force[1] -= rho0*acceCur[1];
        force[2] -= rho0*acceCur[2];

        Cmat *= dvol;
        stre *= dvol;

        for(ii=0;ii<npElem;ii++)
        {
          bb1 = dN_dx[ii];
          bb2 = dN_dy[ii];
          bb3 = dN_dz[ii];
          bb4 = N[ii]*dvol0;
          bb5 = N[ii]*dvol0;

          fact1 = (dNc_dx(ii) - bb1)/3.0;
          fact2 = (dNc_dy(ii) - bb2)/3.0;
          fact3 = (dNc_dz(ii) - bb3)/3.0;

          //fact1 = fact2 = fact3 = 0.0;

          Bbar1(0,0) = bb1+fact1;
          Bbar1(0,1) = 0.0;
          Bbar1(0,2) = 0.0;
          Bbar1(0,3) = bb2;
          Bbar1(0,4) = fact1;
          Bbar1(0,5) = 0.0;
          Bbar1(0,6) = bb3;
          Bbar1(0,7) = 0.0;
          Bbar1(0,8) = fact1;

          Bbar1(1,0) = fact2;
          Bbar1(1,1) = bb1;
          Bbar1(1,2) = 0.0;
          Bbar1(1,3) = 0.0;
          Bbar1(1,4) = bb2+fact2;
          Bbar1(1,5) = 0.0;
          Bbar1(1,6) = 0.0;
          Bbar1(1,7) = bb3;
          Bbar1(1,8) = fact2;

          Bbar1(2,0) = fact3;
          Bbar1(2,1) = 0.0;
          Bbar1(2,2) = bb1;
          Bbar1(2,3) = 0.0;
          Bbar1(2,4) = fact3;
          Bbar1(2,5) = bb2;
          Bbar1(2,6) = 0.0;
          Bbar1(2,7) = 0.0;
          Bbar1(2,8) = bb3+fact3;

          bc = Bbar1*Cmat;

          TI   = 3*ii;
          TIp1 = TI+1;
          TIp2 = TI+2;

          Gxs = Bbar1*stre;

          Flocal(TI)   += (bb5*force[0] - Gxs[0]);
          Flocal(TIp1) += (bb5*force[1] - Gxs[1]);
          Flocal(TIp2) += (bb5*force[2] - Gxs[2]);
          
          for(jj=0; jj<npElem; jj++)
          {
            cc1 = dN_dx[jj];
            cc2 = dN_dy[jj];
            cc3 = dN_dz[jj];
            cc4 = N[jj];

            fact1 = (dNc_dx(jj) - cc1)/3.0;
            fact2 = (dNc_dy(jj) - cc2)/3.0;
            fact3 = (dNc_dz(jj) - cc3)/3.0;

            //fact1 = fact2 = fact3 = 0.0;

            Bbar2(0,0) = cc1+fact1;
            Bbar2(1,0) = 0.0;
            Bbar2(2,0) = 0.0;
            Bbar2(3,0) = cc2;
            Bbar2(4,0) = fact1;
            Bbar2(5,0) = 0.0;
            Bbar2(6,0) = cc3;
            Bbar2(7,0) = 0.0;
            Bbar2(8,0) = fact1;

            Bbar2(0,1) = fact2;
            Bbar2(1,1) = cc1;
            Bbar2(2,1) = 0.0;
            Bbar2(3,1) = 0.0;
            Bbar2(4,1) = cc2+fact2;
            Bbar2(5,1) = 0.0;
            Bbar2(6,1) = 0.0;
            Bbar2(7,1) = cc3;
            Bbar2(8,1) = fact2;

            Bbar2(0,2) = fact3;
            Bbar2(1,2) = 0.0;
            Bbar2(2,2) = cc1;
            Bbar2(3,2) = 0.0;
            Bbar2(4,2) = fact3;
            Bbar2(5,2) = cc2;
            Bbar2(6,2) = 0.0;
            Bbar2(7,2) = 0.0;
            Bbar2(8,2) = cc3+fact3;


            TJ   = 3*jj;
            TJp1 = TJ+1;
            TJp2 = TJ+2;

            fact  = bb4*acceFact1*cc4*rho0;

            Klocal(TI,   TJ)    +=  fact;
            Klocal(TIp1, TJp1)  +=  fact;
            Klocal(TIp2, TJp2)  +=  fact;

            Ktemp = (af*bc)*Bbar2;

            Klocal(TI,   TJ)    += Ktemp(0,0) ;
            Klocal(TI,   TJp1)  += Ktemp(0,1) ;
            Klocal(TI,   TJp2)  += Ktemp(0,2) ;

            Klocal(TIp1, TJ)    += Ktemp(1,0) ;
            Klocal(TIp1, TJp1)  += Ktemp(1,1) ;
            Klocal(TIp1, TJp2)  += Ktemp(1,2) ;

            Klocal(TIp2, TJ)    += Ktemp(2,0) ;
            Klocal(TIp2, TJp1)  += Ktemp(2,1) ;
            Klocal(TIp2, TJp2)  += Ktemp(2,2) ;
          }
        }
    } //gp

    //if(elenum == 0)
      //printMatrix(Klocal);  printf("\n\n\n");  printVector(Flocal);

    return 0;
}




/*
int Elem_Solid_3D_HEXA8BbarFbar::calcStiffnessAndResidualFS(MatrixXd& Klocal, VectorXd& Flocal)
{
    int   err,  isw,  count,  count1, index, ii, jj, kk, ll, mm, gp, TI, TIp1, TIp2, TJ, TJp1, TJp2;

    VectorXd  N(npElem), dN_dx(npElem), dN_dy(npElem), dN_dz(npElem);

    double  detF, detFc, trF, trFc, fact, fact1, dvol, dvol0, Jac, volstrain, r2d3 = 2.0/3.0;
    double  bb1, bb2, bb3, bb4, bb5, cc1, cc2, cc3, cc4, cc5;
    double  param[3], bforce[3], force[3], tarr[6];
    double  veloCur[3], acceCur[3], sig[3];
    MatrixXd  F(3,3), Fn(3,3), Fc(3,3), Cmat(9,9), Bmat(3,9), Gc(3,9);
    Bmat.setZero();
    Cmat.setZero();
    Gc.setZero();
    F.setZero();
    Fn.setZero();
    Fc.setZero();
    VectorXd  stre(9), qVec(9);
    stre.setZero();

    bool  finite = MatlData->isFiniteStrain();

    double rho0  = MatlData->getDensity() ;
    double rho   = rho0;

    //bforce[0]   = elmDat[6]*timeFunction[0]->getFactor() ;
    //bforce[1]   = elmDat[7]*timeFunction[0]->getFactor() ;
    //bforce[2]   = elmDat[8]*timeFunction[0]->getFactor() ;

    double af   = SolnData->td(2);
    double acceFact1 = SolnData->td(5);
    double acceFact2 = acceFact1;
    double FiniteFact = (finite == true) ? 1.0 : 0.0;

    double xNode[npElem], yNode[npElem], zNode[npElem], xx, yy,  zz;
    for(ii=0;ii<npElem;ii++)
    {
      xNode[ii] = GeomData->NodePosOrig[nodeNums[ii]][0];
      yNode[ii] = GeomData->NodePosOrig[nodeNums[ii]][1];
      zNode[ii] = GeomData->NodePosOrig[nodeNums[ii]][2];
    }


    int nGP=1;
    vector<double>  gausspoints1, gausspoints2, gausspoints3, gaussweights;

    // Calculate values at the centroid
    if(npElem == 6)
    {
      getGaussPointsWedge(nGP, gausspoints1, gausspoints2, gausspoints3, gaussweights);
    }
    else if(npElem == 8)
    {
      getGaussPointsHexa(nGP, gausspoints1, gausspoints2, gausspoints3, gaussweights);
    }
    else if(npElem == 10)
    {
      getGaussPointsTetra(nGP, gausspoints1, gausspoints2, gausspoints3, gaussweights);
    }
    else if(npElem == 18)
    {
      getGaussPointsWedge(nGP, gausspoints1, gausspoints2, gausspoints3, gaussweights);
    }
    else if(npElem == 27)
    {
      getGaussPointsHexa(nGP, gausspoints1, gausspoints2, gausspoints3, gaussweights);
    }


    param[0] = gausspoints1[gp];
    param[1] = gausspoints2[gp];
    param[2] = gausspoints3[gp];

    GeomData->computeBasisFunctions3D(CONFIG_ORIGINAL, ELEM_SHAPE, degree, param, nodeNums, &N(0), &dNc_dx(0), &dNc_dy(0), &dNc_dz(0), Jac);

    computeDefGradCur(dNc_dx, dNc_dy, dNc_dz, Fc);

    trFc  = Fc.trace();
    detFc = Fc.determinant();

    // basis functions wrt to the current configuration
    GeomData->computeBasisFunctions3D(CONFIG_DEFORMED, ELEM_SHAPE, degree, param, nodeNums, &N(0), &dNc_dx(0), &dNc_dy(0), &dNc_dz(0), Jac);



    if(npElem == 6)
    {
      nGP = 2;
      getGaussPointsWedge(nGP, gausspoints1, gausspoints2, gausspoints3, gaussweights);
    }
    else if(npElem == 8)
    {
      nGP = 8;
      getGaussPointsHexa(nGP, gausspoints1, gausspoints2, gausspoints3, gaussweights);
    }
    else if(npElem == 10)
    {
      nGP = 4;
      getGaussPointsTetra(nGP, gausspoints1, gausspoints2, gausspoints3, gaussweights);
    }
    else if(npElem == 18)
    {
      nGP = 9;
      getGaussPointsWedge(nGP, gausspoints1, gausspoints2, gausspoints3, gaussweights);
    }
    else if(npElem == 27)
    {
      nGP = 27;
      getGaussPointsHexa(nGP, gausspoints1, gausspoints2, gausspoints3, gaussweights);
    }



    if(Klocal.rows() != nsize)
    {
      Klocal.resize(nsize, nsize);
      Flocal.resize(nsize);
    }
    Klocal.setZero();
    Flocal.setZero();


    for(gp=0; gp<nGP; gp++)
    {
        param[0] = gausspoints1[gp];
        param[1] = gausspoints2[gp];
        param[2] = gausspoints3[gp];

        GeomData->computeBasisFunctions3D(CONFIG_ORIGINAL, ELEM_SHAPE, degree, param, nodeNums, &N(0), &dN_dx(0), &dN_dy(0), &dN_dz(0), Jac);

        dvol0 = gaussweights[gp]*Jac;
        dvol  = dvol0;

        computeDefGradCur(dN_dx, dN_dy, dN_dz, F);

        volstrain = F.trace() - 3.0;

        detF = F.determinant();

        if(detF < 0.0)
        {
          throw runtime_error("Negative Jacobian in the element");
        }

        //acceCur[0] = computeValueDotDotCur(0, N);
        //acceCur[1] = computeValueDotDotCur(1, N);
        //acceCur[2] = computeValueDotDotCur(2, N);

        if(finite)
        {
          GeomData->computeBasisFunctions3D(CONFIG_DEFORMED, ELEM_SHAPE, degree, param, nodeNums, &N(0), &dN_dx(0), &dN_dy(0), &dN_dz(0), Jac);
          dvol = gaussweights[gp]*Jac;
        }

        //  modify F
        fact = pow(detFc/detF, 1.0/3.0);
        F *= fact ;

        MatlData->computeStressAndTangent(true, sss, F, F, pres, stre, Cmat, ivar, gp, myTime.dt);

        // q tensor for the F-bar formulation
        qVec = -r2d3*stre;
        for(ii=0; ii<9; ii++)
        {
            qVec(ii) += ( Cmat(ii,0)+Cmat(ii,4)+Cmat(ii,8)+stre(ii) )/3.0;
        }


        // Calculate Stiffness and Residual
        //==============================================

        // body force and acceleration terms
        force[0] = rho0*bforce[0];
        force[1] = rho0*bforce[1];
        force[2] = rho0*bforce[2];

        force[0] -= rho0*acceCur[0];
        force[1] -= rho0*acceCur[1];
        force[2] -= rho0*acceCur[2];

        for(ii=0; ii<npElem; ii++)
        {
          bb1 = dvol*dN_dx[ii];
          bb2 = dvol*dN_dy[ii];
          bb3 = dvol*dN_dz[ii];
          bb4 = dvol*N[ii];
          bb5 = dvol0*N[ii];

          TI   = 3*ii;
          TIp1 = TI+1;
          TIp2 = TI+2;

          //  Klocal matrix
          for(jj=0; jj<9; jj++)
          {
            Gc(0,jj) = bb1*Cmat(0,jj) + bb2*Cmat(3,jj) + bb3*Cmat(6,jj);
            Gc(1,jj) = bb1*Cmat(1,jj) + bb2*Cmat(4,jj) + bb3*Cmat(7,jj);
            Gc(2,jj) = bb1*Cmat(2,jj) + bb2*Cmat(5,jj) + bb3*Cmat(8,jj);
          }

          sig[0] = bb1*stre[0] + bb2*stre[3] + bb3*stre[6];
          sig[1] = bb1*stre[1] + bb2*stre[4] + bb3*stre[7];
          sig[2] = bb1*stre[2] + bb2*stre[5] + bb3*stre[8];

          Flocal(TI)   += (bb5*force[0] - sig[0]) ;
          Flocal(TIp1) += (bb5*force[1] - sig[1]) ;
          Flocal(TIp2) += (bb5*force[2] - sig[2]) ;

          tarr[0] = bb1*qVec[0] + bb2*qVec[3] + bb3*qVec[6] ;
          tarr[1] = bb1*qVec[1] + bb2*qVec[4] + bb3*qVec[7] ;
          tarr[2] = bb1*qVec[2] + bb2*qVec[5] + bb3*qVec[8] ;

          for(jj=0; jj<npElem; jj++)
          {
            cc1 = dN_dx[jj];
            cc2 = dN_dy[jj];
            cc3 = dN_dz[jj];
            cc4 = N[jj];

            TJ   = 3*jj;
            TJp1 = TJ+1;
            TJp2 = TJ+2;

            // acceleration term
            fact  = bb5*acceFact1*cc4*rho0;

            // material Stiffness contribution
            //fact  += af*(sig[0]*cc1+sig[1]*cc2+sig[2]*cc3)*FiniteFact;

            Klocal(TI,   TJ)    += fact ;
            Klocal(TIp1, TJp1)  += fact ;
            Klocal(TIp2, TJp2)  += fact ;

            Klocal(TI,   TJ)    +=  af*(Gc(0,0)*cc1 + Gc(0,3)*cc2 + Gc(0,6)*cc3) ;
            Klocal(TI,   TJp1)  +=  af*(Gc(0,1)*cc1 + Gc(0,4)*cc2 + Gc(0,7)*cc3) ;
            Klocal(TI,   TJp2)  +=  af*(Gc(0,2)*cc1 + Gc(0,5)*cc2 + Gc(0,8)*cc3) ;

            Klocal(TIp1, TJ)    +=  af*(Gc(1,0)*cc1 + Gc(1,3)*cc2 + Gc(1,6)*cc3) ;
            Klocal(TIp1, TJp1)  +=  af*(Gc(1,1)*cc1 + Gc(1,4)*cc2 + Gc(1,7)*cc3) ;
            Klocal(TIp1, TJp2)  +=  af*(Gc(1,2)*cc1 + Gc(1,5)*cc2 + Gc(1,8)*cc3) ;

            Klocal(TIp2, TJ)    +=  af*(Gc(2,0)*cc1 + Gc(2,3)*cc2 + Gc(2,6)*cc3) ;
            Klocal(TIp2, TJp1)  +=  af*(Gc(2,1)*cc1 + Gc(2,4)*cc2 + Gc(2,7)*cc3) ;
            Klocal(TIp2, TJp2)  +=  af*(Gc(2,2)*cc1 + Gc(2,5)*cc2 + Gc(2,8)*cc3) ;

            //  additional terms due to F-bar formulation
            cc1 = af*(dNc_dx(jj) - cc1);
            cc2 = af*(dNc_dy(jj) - cc2);
            cc3 = af*(dNc_dz(jj) - cc3);

            Klocal(TI,   TJ)    +=  tarr[0]*cc1;
            Klocal(TI,   TJp1)  +=  tarr[0]*cc2;
            Klocal(TI,   TJp2)  +=  tarr[0]*cc3;

            Klocal(TIp1, TJ)    +=  tarr[1]*cc1;
            Klocal(TIp1, TJp1)  +=  tarr[1]*cc2;
            Klocal(TIp1, TJp2)  +=  tarr[1]*cc3;

            Klocal(TIp2, TJ)    +=  tarr[2]*cc1;
            Klocal(TIp2, TJp1)  +=  tarr[2]*cc2;
            Klocal(TIp2, TJp2)  +=  tarr[2]*cc3;
          }
        }
    }//gp

    return 0;
}
*/






int Elem_Solid_3D_HEXA8::calcStiffnessAndResidual_Formulation_Mixed(MatrixXd& Klocal, VectorXd& Flocal)
{
    return 0;
}









int Elem_Solid_3D_HEXA8::calcResidual(VectorXd& Flocal)
{

    return 0;
}





void Elem_Solid_3D_HEXA8::elementContourplot(int vartype, int varindex, int index)
{
   double outval[50];

   switch(vartype)
   {
       case 1:  // plot total strain
       case 2:  // plot elastic strain
       case 3:  // plot plastic strain

                projectStrain(true, vartype, varindex, index, outval);

              break;

       case 4:  // plot stress

                projectStress(true, vartype, varindex, index, outval);

              break;

       case 5:  // plot element internal variables

                projectInternalVariable(true, vartype, varindex, index, outval);

              break;

       default:

              cout  << " Invalid Variable Type to project in 'BernsteinElem2DSolidTria6Node::projectToNodes'" << endl;
              break;
    }

    assert(vals2project.size() >= 1);

    vals2project[0] = 0.0;
    for(int ii=0; ii<nGP; ii++)
      vals2project[0] += outval[ii];

    vals2project[0] /= nGP;

    return;
}


void Elem_Solid_3D_HEXA8::projectToNodes(bool extrapolateFlag, int vartype, int varindex, int index)
{
   double outval[50];

   switch(vartype)
   {
       case 1:  // plot total strain
       case 2:  // plot elastic strain
       case 3:  // plot plastic strain

                projectStrain(extrapolateFlag, vartype, varindex, index, outval);

              break;

       case 4:  // plot stress

                projectStress(extrapolateFlag, vartype, varindex, index, outval);

              break;

       case 5:  // plot element internal variables

                projectInternalVariable(extrapolateFlag, vartype, varindex, index, outval);

              break;

       default:

              cout  << "           Invalid Variable Type to project " << endl;
              break;
    }


    assert(vals2project.size() == npElem);
    stressrecovery_extrapolate_Hexahedron(npElem, outval, &vals2project[0]);

  return;
}


void Elem_Solid_3D_HEXA8::projectStress(bool extrapolateFlag, int vartype, int varindex, int index, double* outval)
{


    return;
}



void Elem_Solid_3D_HEXA8::projectStrain(bool extrapolateFlag, int vartype, int varindex, int index, double* outval)
{
  return;
}



void Elem_Solid_3D_HEXA8::projectInternalVariable(bool extrapolateFlag, int vartype, int varindex, int index, double* outval)
{
  return;
}




void  Elem_Solid_3D_HEXA8::computeEnergy(int index1, int index2, VectorXd& energy)
{
    // compute total energy for structural dynamic problems
    ///////////////////////////////////////////////////////////

    return;
}


